home *** CD-ROM | disk | FTP | other *** search
/ Graphics Plus / Graphics Plus.iso / general / modelers / geomview / source.lha / Geomview / src / lib / gprim / bezier / bezdice.c < prev    next >
Encoding:
C/C++ Source or Header  |  1993-07-20  |  8.3 KB  |  280 lines
  1. /* Copyright (c) 1992 The Geometry Center; University of Minnesota
  2.    1300 South Second Street;  Minneapolis, MN  55454, USA;
  3.    
  4. This file is part of geomview/OOGL. geomview/OOGL is free software;
  5. you can redistribute it and/or modify it only under the terms given in
  6. the file COPYING, which you should have received along with this file.
  7. This and other related software may be obtained via anonymous ftp from
  8. geom.umn.edu; email: software@geom.umn.edu. */
  9.  
  10. /* Authors: Silvio Levy, Charlie Gunn, Stuart Levy, Tamara Munzner, 
  11.    Mark Phillips */
  12.  
  13. #include "math.h"
  14. #include "bezierP.h"
  15.  
  16. Bezier *
  17. BezierDice( Bezier *bezier, int nu, int nv )
  18. {
  19.     if(nu < 2) nu = BEZ_DEFAULT_MESH_SIZE;
  20.     if(nv < 2) nv = BEZ_DEFAULT_MESH_SIZE;
  21.     if(nu != bezier->nu || nv != bezier->nv) {
  22.     bezier->flag |= BEZ_REMESH;
  23.     bezier->nu = nu;
  24.     bezier->nv = nv;
  25.     }
  26.     return bezier;
  27. }
  28.  
  29. /* this gets called when it's time to remesh, so the current mesh is destroyed
  30.    and a new one is allocated */
  31. Mesh *
  32. BezierReDice(bezier)
  33.     Bezier *bezier;
  34. {
  35.     int      u, v, d;
  36.     int      dimn, nu, nv, degree_u, degree_v;
  37.     float    *tmpp0, *tmpp1, *tmpp2, *p;
  38.     float    *tmpdv, *tmpdv0, *tmpdv1, *tmpdv2, *dv;
  39.     float    *tmpdu, *tmpdu0, *tmpdu1, *tmpdu2, *du;
  40.     float norm;
  41.     HPoint3   *bmp;
  42.     Point3 *bmn;
  43.     Mesh m;
  44.  
  45.     dimn = bezier->dimn;
  46.     degree_u = bezier->degree_u;
  47.     degree_v = bezier->degree_v;
  48.  
  49.     if (bezier->nu <= 1)         bezier->nu = BEZ_DEFAULT_MESH_SIZE;
  50.     if (bezier->nv <= 1)         bezier->nv = BEZ_DEFAULT_MESH_SIZE;
  51.     nu = bezier->nu;
  52.     nv = bezier->nv;
  53.  
  54.     m.nu = nu;
  55.     m.nv = nv;
  56.     m.flag = 0;
  57.     m.geomflags = 0;
  58.     m.flag |= MESH_N;
  59.     if (bezier->flag & BEZ_C)    m.flag |= MESH_C;
  60.     m.p = NULL;
  61.     m.c = NULL;
  62.     m.n = NULL;
  63.  
  64.     m.p = OOGLNewNE(HPoint3, nu*nv, "BezierReDice: mesh points");
  65.     m.n = OOGLNewNE(Point3, nu*nv, "BezierReDice: mesh normals");
  66.     if(bezier->flag & BEZ_C)
  67.       m.c = OOGLNewNE(ColorA, nu*nv,  "BezierReDice: mesh colors");
  68.  
  69.  
  70.     /* compute first pass: interpolate mesh points in u direction */
  71.     tmpdv = OOGLNewNE(float, dimn * degree_v * (degree_u+1),
  72.               "BezierReDice: tempdv");
  73. /*
  74. printf("this is dv\n");
  75. */
  76.     for (v=0; v<degree_v; ++v)
  77.       for (u=0; u<degree_u+1; ++u) {
  78.         for (d=0; d<dimn; ++d) {
  79.           *(tmpdv+dimn*(v*(degree_u+1)+u)+d)=
  80.             *(bezier->CtrlPnts + dimn * ((v+1)*(degree_u+1)+u)+d)-
  81.             *(bezier->CtrlPnts + dimn * (v*(degree_u+1)+u)+d);
  82. /*
  83.           printf("%f ", *(tmpdv+dimn*(v*(degree_u+1)+u)+d));
  84. */
  85.         }
  86. /*
  87.       printf("\n");
  88. */
  89.       }
  90.        
  91.     tmpdu = OOGLNewNE(float,dimn*(degree_v+1)*degree_u,"BezierReDice:tmpdu");
  92. /*
  93. printf("this is du\n");
  94. */
  95.     for (v=0; v<degree_v+1; ++v)
  96.       for (u=0; u<degree_u; ++u) {
  97.         for (d=0; d<dimn; ++d) {
  98.           *(tmpdu+dimn*(v*degree_u+u)+d)=
  99.             *(bezier->CtrlPnts + dimn * (v*(degree_u+1)+u+1)+d)-
  100.             *(bezier->CtrlPnts + dimn * (v*(degree_u+1)+u)+d);
  101. /*
  102.           printf("%f ", *(tmpdv+dimn*(v*degree_u+u)+d));
  103. */
  104.         }
  105. /*
  106.       printf("\n");
  107. */
  108.       }
  109.     
  110. /*
  111. printf("this is tmpp0\n");
  112. */
  113.     tmpp0 = OOGLNewNE(float, dimn * (degree_v+1) * nu, "BezierReDice: tmpp0");
  114.     for (v=0; v<degree_v+1; ++v)
  115.       bezier_interp(bezier->CtrlPnts + v * dimn * (degree_u+1),
  116.                     tmpp0 + v * dimn * nu, degree_u, nu, dimn);
  117.     tmpdv0 = OOGLNewNE(float, dimn * degree_v * nu, "BezierReDice: tmpdv0");
  118. /*
  119. printf("this is tmpdv0\n");
  120. */
  121.     for (v=0; v<degree_v; ++v)
  122.       bezier_interp(tmpdv + v * dimn * (degree_u+1),
  123.                     tmpdv0 + v * dimn * nu, degree_u, nu, dimn);
  124. /*
  125. printf("this is tmpdu0\n");
  126. */
  127.     tmpdu0 = OOGLNewNE(float, dimn * (degree_v+1) * nu,"BezierReDice: tmpdu0");
  128.     for (v=0; v<degree_v+1; ++v)
  129.       bezier_interp(tmpdu + v * dimn * degree_u,
  130.                     tmpdu0 + v * dimn * nu, degree_u-1, nu, dimn);
  131.  
  132.     /* now compute second pass, filling in columns */
  133.     tmpp1 = OOGLNewNE(float, dimn * (degree_v+1), "BezierReDice: tmpp1");
  134.     tmpdv1 = OOGLNewNE(float, dimn * degree_v, "BezierReDice: tmpdv1");
  135.     tmpdu1 = OOGLNewNE(float, dimn * (degree_v+1), "BezierReDice: tmpdu1");
  136.     tmpp2 = OOGLNewNE(float, dimn * nv, "BezierReDice: tmpp2");
  137.     tmpdv2 = OOGLNewNE(float, dimn * nv, "BezierReDice: tmpdv2");
  138.     tmpdu2 = OOGLNewNE(float, dimn * nv, "BezierReDice: tmpdu2");
  139.     bmp = m.p;
  140.     bmn = m.n;
  141.     for (u=0; u<nu; ++u) {
  142.       for (v=0; v<degree_v+1; v++)
  143.         bcopy(tmpp0 + (u + v*nu)*dimn, tmpp1+v*dimn, dimn*sizeof(float));
  144.       for (v=0; v<degree_v; v++)
  145.         bcopy(tmpdv0 + (u + v*nu)*dimn, tmpdv1+v*dimn, dimn*sizeof(float));
  146.       for (v=0; v<degree_v+1; v++)
  147.         bcopy(tmpdu0 + (u + v*nu)*dimn, tmpdu1+v*dimn, dimn*sizeof(float));
  148.  
  149.       bezier_interp(tmpp1, tmpp2, degree_v, nv, dimn);
  150.       bezier_interp(tmpdv1, tmpdv2, degree_v-1, nv, dimn);
  151.       bezier_interp(tmpdu1, tmpdu2, degree_v, nv, dimn);
  152.       for (v=0, p=tmpp2, dv=tmpdv2, du=tmpdu2; 
  153.            v<nv; ++v, p+=dimn, dv+=dimn, du+=dimn) {
  154.         if(dimn == 4) {
  155.           dv[0] = dv[0]*p[3]-dv[3]*p[0];
  156.           dv[1] = dv[1]*p[3]-dv[3]*p[1];
  157.           dv[2] = dv[2]*p[3]-dv[3]*p[2];
  158.           du[0] = du[0]*p[3]-du[3]*p[0];
  159.           du[1] = du[1]*p[3]-du[3]*p[1];
  160.           du[2] = du[2]*p[3]-du[3]*p[2];
  161.           p[0] /= p[3];
  162.           p[1] /= p[3];
  163.           p[2] /= p[3];
  164.         }
  165.         *bmp = *(HPoint3 *)p;   bmp->w = 1.; bmp++;
  166.         bmn->x = -du[1]*dv[2]+du[2]*dv[1];
  167.         bmn->y = -du[2]*dv[0]+du[0]*dv[2];
  168.         bmn->z = -du[0]*dv[1]+du[1]*dv[0];
  169.         norm=bmn->x *bmn->x +bmn->y *bmn->y +bmn->z *bmn->z;
  170.         if (norm==0.0)
  171.           {bmn->x=1.; norm=1.;}
  172.         norm=sqrt(norm);
  173.         bmn->x /= norm;
  174.         bmn->y /= norm;
  175.         bmn->z /= norm;
  176.         bmn++;
  177.       }
  178.     }
  179.     GeomFree(tmpp0);
  180.     GeomFree(tmpp1);
  181.     GeomFree(tmpp2);
  182.     GeomFree(tmpdu);
  183.     GeomFree(tmpdu0);
  184.     GeomFree(tmpdu1);
  185.     GeomFree(tmpdu2);
  186.     GeomFree(tmpdv);
  187.     GeomFree(tmpdv0);
  188.     GeomFree(tmpdv1);
  189.     GeomFree(tmpdv2);
  190.     bezier->flag &= ~BEZ_REMESH;        /* turn off this bit */
  191.  
  192.     if(bezier->flag & BEZ_C) {
  193.         float fu, unu, fv, unv;
  194.         ColorA u0, u1;
  195.         ColorA *cp;
  196.  
  197. #define INTC(c0, c1, t, unt, dest)        /* Interpolate color */ \
  198.                 dest.r = c0.r * unt + c1.r * t; \
  199.                 dest.g = c0.g * unt + c1.g * t; \
  200.                 dest.b = c0.b * unt + c1.b * t; \
  201.                 dest.a = c0.a * unt + c1.a * t;
  202.  
  203.         cp = m.c;
  204.         for(v = 0; v < nv; v++) {
  205.            fv = (float) v / (nv - 1);
  206.            unv = 1 - fv;
  207.            INTC(bezier->c[0], bezier->c[1], fv, unv, u0);
  208.            INTC(bezier->c[2], bezier->c[3], fv, unv, u1);
  209.            for(u = 0; u < nu; u++) {
  210.                 fu = (float) u / (nu - 1);
  211.                 unu = 1 - fu;
  212.                 INTC(u0, u1, fu, unu, (*cp));
  213.                 cp++;
  214.            }
  215.         }
  216.     }
  217.     GeomDelete((Geom *)bezier->mesh);
  218.     if((bezier->mesh = (Mesh *) GeomCreate("mesh",
  219.         CR_NOCOPY,
  220. /*     For now, assume these are rational beziers hence belong in 3-space
  221.         CR_4D, (dimn == 4) ? 1 : 0,
  222. */
  223.         CR_FLAG,m.flag,
  224.         CR_NU,m.nu,
  225.         CR_NV,m.nv, 
  226.         CR_POINT4,m.p,
  227.         CR_NORMAL,m.n,
  228.         CR_COLOR,m.c,
  229.         CR_END)) == NULL) {
  230.     OOGLError(1, "BezierReDice: can't create Mesh");
  231.     return NULL;
  232.     }
  233.  
  234.     return(bezier->mesh);
  235. }
  236.  
  237.  
  238. /* Takes a list of control points describing a spline */
  239. /* and produces values for the spline */
  240. /* at  n  equally spaced points.   */
  241.  
  242. bezier_interp(in,out,deg,n,dimn)
  243. float *in; /* input array of control points */
  244. float *out; /* output array */
  245. int deg; /* degree of spline */
  246. int n; /* number of points to interpolate */
  247. int dimn; /* dimension of range */
  248. {
  249.   int    m, offset, k, j;
  250.   float  p[(MAX_BEZ_DEGREE+1)*MAX_BEZ_DIMN];
  251.   float  t;
  252.  
  253. /*
  254. for (j=0; j<deg+1; ++j) 
  255.   printf("%f %f %f    ",*(in+j*dimn),*(in+j*dimn+1),*(in+j*dimn+2));
  256.   printf("\n");
  257. */
  258.  
  259.   for (j=0; j<n; ++j) {
  260.     t = ((float) j) / (n - 1);
  261.     bcopy(in, p, dimn*sizeof(float)*(deg+1));
  262.     for (k=0; k<deg; ++k)
  263.       for (offset=0, m=0; m<deg; ++m, offset += dimn) {
  264.         p[offset] += t * (p[offset+dimn] - p[offset]);
  265.         p[offset+1] += t * (p[offset+dimn+1] - p[offset+1]);
  266.         p[offset+2] += t * (p[offset+dimn+2] - p[offset+2]);
  267.         if(dimn == 4)
  268.            p[offset+3] += t * (p[offset+dimn+3] - p[offset+3]);
  269.     }
  270.     bcopy(p, out, sizeof(float)*dimn);
  271.     out+=dimn;
  272.   }
  273. /*
  274. out -= dimn*n;
  275. for (j=0; j<n; ++j) 
  276.   printf("%f %f %f    ",*(out+j*dimn),*(out+j*dimn+1),*(out+j*dimn+2));
  277.   printf("\n");
  278. */
  279. }
  280.